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^ : Abstract 

oo 

A transfer matrix scaling technique is developed for randomly diluted sys- 



tems, and applied to the site-diluted Ising model on a square lattice in two 
dimensions. For each allowed disorder configuration between two adjacent 



C^^ \ columns, the contribution of the respective transfer matrix to the decay of 

1^ I correlations is considered only as far as the ratio of its two largest eigenval- 

f^ ■ ues, allowing an economical calculation of a configuration-averaged correlation 

^J . length. Standard phenomenological-renormalisation procedures are then used 

to analyse aspects of the phase boundary which are difficult to assess accu- 
C^ \ rately by alternative methods. For magnetic site concentration p close to p, 



C) 



the extent of exponential behaviour of the Tc x p curve is clearly seen for 



'^ ■ over two decades of variation of p — Pc- Close to the pure-system limit, the 

^ ■ exactly-known reduced slope is reproduced to a very good approximation, 

Q . though with non-monotonic convergence. The averaged correlation lengths 

are inserted into the exponent-amplitude relationship predicted by conformal 
invariance to hold at criticality. The resulting exponent r] remains near the 
;h ■ pure value (1/4) for all intermediate concentrations until it crosses over to the 

percolation value at the threshold. 
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The transfer matrix is well known to provide exact solutions for a number of low- 
dimensional pure systems such as spin models etc [|[]. Applied to finite- width strips and 
combined with finite size scaling 0, it has also proved to be extremely powerful and accu- 
rate for non-solvable pure cases, ranging from magnetic systems |0 to walks [Q and quantum 
Hamiltonians |Q . The application of such "strip-scaling" techniques to random systems has 
however been extremely limited, despite the enormous interest in such problems as the spin 
glass, random field and dilute magnets, ceramic superconductors etc. This is because the 
form so far utilized [§-|^ applies the strip scaling to particular realisations of the random 
system which need not be representative unless appropriate averaging (or self-averaging) is 
employed, resulting in very large scale computing on extremely long strips and/or many 
realisations. 

Here, we provide a strip scaling approach for random systems, in which the disorder 
averaging is carried out as one proceeds along the strip, and which therefore does not have the 
deficiencies noted above. Results reported here extend, and give details of, those contained 
in a previous Rapid Communication |T^. Our scheme does rely on certain assumptions 
regarding the dominant contributions to the decay of correlations, which must be spelt out 
clearly and supported by numerical evidence whenever possible. This is one of our purposes 
in what follows. 

The main numerical application of the technique here is to the two-dimensional site- 
diluted Ising model ||ri|. The transfer- matrix descriptions of the special limiting cases of 
percolation |jl2[ and pure Ising system |T^ are well-understood, and are reproduced by the 



present approach, as we shall see below. On the other hand, several aspects of the behaviour 
at intermediate dilution are either not known to the same degree of accuracy, or have not 
been successfully accounted for via transfer-matrix methods. Of the former, we shall focus 
on the extent of exponential behaviour of the T^ x p curve for magnetic site concentration 
p close to pc, and on the variation of the decay-of-correlations exponent rj along the critical 
line; as for the latter, we give results for the reduced slope of the phase boundary close to 



the pure-system limit, to be compared with an exact solution by Thorpe and McGurn |]I4 

In the transfer-matrix scaling formulation, the crucial quantity to be considered is the 
rate of decay of correlations. On a strip, the correlation functions always (i.e. at finite 
temperatures) decay exponentially for suitably large distances: 

< aoCTR >~ exp(-i?/0 , (1) 

with the correlation length C, being given, for a pure system, by 

r' = -lnA2/Ai, (2) 



where Ai, A2 are the two largest eigenvalues of the transfer matrix [15|. To obtain this pure 
system result, the effect of repeated applications of the (same) transfer matrix is investigated, 
and only leading terms are retained, again in the long-distance limit. In order to have a 
quantitative understanding of this limit, it is useful to calculate the correlation functions 
explicitly for the pure Ising model on strips, and compare their actual rate of decay to that 
predicted by Eqs. [| and ||. This is illustrated in Fig. |I|, where two typical examples are 
shown. One can see that, for widely differing temperatures and strip widths, the influence 
of subdominant terms on the correlation decay is always very small and becomes negligible 
after no more than two lattice spacings. We shall invoke this in what follows. 



For a random system, the decay of thermal correlations for a given realisation of random- 
ness is still given through the iterative product of transfer matrices from one spin column 
to the next; however, some form of configurational averaging is now necessary. In previous 
work [§-§] this was done by generating strips of length iV >> 1 and carrying out the ac- 
tual transfer-matrix products, so the end result would presumably reflect the properties of 
a representative sample. While the average free energy per site is unambiguously related 
to the dominant Lyapunov exponent of the transfer- matrix product |]TB[, the relationship 



of correlation functions to Lyapunov exponents is much more involved, because the prob- 
ability distribution of the correlation functions themselves turns out to be rather complex 
|T^,|T8|. It can be shown that the most probable correlation decay is given by the difference 
between the two largest Lyapunov exponents [0, which constitutes a straightforward gen- 
eralisation of Eq. 0. This may not be the same as the average correlation decay, which is 
the experimentally measurable quantity; the distinction between the two becomes especially 
important when the (effective) interactions fluctuate in sign, as is the case of spin glasses 



or when random fields are present |[T7|JT8[1 . However, when all interactions are ferromagnetic 
numerical evidence |]T9[ indicates that, at least for two-dimensional Ising systems, the two 
quantities yield essentially the same results apart from logarithmic corrections. In what 
follows, we shall assume this to be the case also when there are absent sites, all non-zero 
interactions being positive. 

The implementation of the schemes just described runs into the usual problems of judg- 
ing when the strip is sufficiently long, so that the sample can be taken as representative. 
For diluted systems the additional difficulty would arise of how to avoid the effects of dis- 



connections. In early attempts |^0|] this was done by introducing a fictitious "weak " bond 
(of strength 10~^ that of "normal" ones) , linking pairs of sites of which at least one was 
absent . Numerical results were, however, unsatisfactory. 

Here, we perform the averages in a different way. Our starting point is the transfer-matrix 
formulation of the percolation problem |T^. In this case, which is the zero-temperature limit 



of diluted magnets, the (geometric) correlation length is correctly given by the largest eigen- 
value of a transfer matrix, whose elements are related to the probability that two adjacent 
site columns (each with occupied and vacant sites) are linked to each other and to the 
origin. Note that disconnections are altogether avoided from the start, as only column con- 
figurations where at least one occupied site is part of the infinite cluster |]I2| are considered. 
Unconnected configurations can be discarded, as they just account for conservation of total 
probability. Were they to be included, the matrix would have 1 as its largest eigenvalue, 
with the second largest eigenvalue giving the correlation length, in exact correspondence 
with Eq. 1^ 1^ Care is thus taken of the correlations between frozen spins in the dilute spin 
system at absolute zero; as the temperature is raised from zero, each of the possible links 
represented by the non-zero elements of the geometric transfer matrix is weakened by ther- 
mal fiuctuations. We propose to take these into account still within the framework of a single 
matrix, so the exact column-by-column character of the average of disorder configurations 
will be preserved. 

Note however that any approach which directly combines a thermal and disorder av- 
eraging (in e.g. multiplying an extended transfer matrix) will be physically incorrect for 
quenched disorder problems. 

So, we must suitably modify each non-zero element of the original matrix; obviously we 



must look at the properties of the spin transfer matrix between the corresponding site column 
states. As the physical property under investigation is the rate of decay of correlations, we 
draw on an analogy to the leading contribution towards this in periodic systems: we choose 
to take the ratio of the two largest eigenvalues of the spin transfer matrix and multiply the 
corresponding geometric transfer matrix element by it. If z and j are two site column states 
connected to each other and to the origin, with respective probabilities Pi and Pj, and the 
spin transfer matrix T*-' between these columns has A*/ and X2 ^-^ its largest eigenvalues, the 
matrix element of the "thermal-geometric" transfer matrix of our formulation is then : 



%, = ^P,P,{Xi/XI) . (3) 

The averaged correlation length in this approximation is given by the largest eigenvalue, Ai, 
of T via 

(r')a.e = -lnAi . (4) 

The following comments are in order: 

(a) As the temperature T — > 0, the two largest eigenvalues of all thermal transfer matrices 
become degenerate, and Eq. R shows that T reduces to the geometric transfer matrix of 



Ref. |I2|, as it should. Thus, the present calculational scheme is asymptotically correct 
in the low-temperature regime. For other disordered systems this must be true as well, 
provided that one can start from a suitable transfer-matrix description of the ground-state 
correlations. 

(b) As the concentration p of magnetic sites approaches 1, the only remaining non-zero 
element of T is along the diagonal, connecting two fully occupied columns; denoting by \p^^^ 
and X^"^^ the largest eigenvalues of the corresponding thermal transfer matrix, 

Ai = ArvAr', p = 1 (5) 

and the pure system limit is correctly obtained. 

(c) As befits quenched problems, disorder and thermal aspects are not being averaged 
together. We represent the effect of thermal fluctuations on each given geometric configu- 
ration by the ratio of the respective largest thermal eigenvalues, and the disorder average 
is performed at a later step, in obtaining the largest eigenvalue of T. Thus we make the 
analogue of the configurational average of the factor e~^/^ in the correlation function (related 
to its decay between two adjacent columns) and it is indeed the correlation function which 
is self-averaging. 

(d) The procedure outlined contains the approximation that the contribution given by 
each thermal transfer matrix is truncated and replaced by that of its two largest eigenvalues. 
Once this has been done, we follow lines analogous to those of Refs. |^J^,|] for the calculation 
of an averaged correlation length in a random system. We are thus (approximately) aver- 
aging the two largest Lyapunov exponents, and not the correlation functions. If all thermal 
transfer matrices, corresponding to distinct disorder configurations, commuted with each 
other (which is not the case) our ansatz would be identical to that of Refs. 0,^,||; note that 
the probabilistic factors included in our matrix elements would be mirrored in the relative 
frequency of configurations generated by those methods. On the other hand, the influence 
of higher-ranking eigenvalues is expected to die out asymptotically, thus minimising the 



truncation effects mentioned above. Having in mind the corresponding results for the pure 
Ising case, portrayed in Fig |l| above, this does not seem too drastic an assumption. 

The results for the critical curve are now presented. 

Along the critical line T^ = Tc{p), the correlation length ^ diverges in the infinite system 
limit. In phenomenological scaling between strip widths L and L' this is where |^ ^l/L = 
C,L'/L'. This condition would normally give a fixed point in a one-parameter space. Here, 
where two variables T, p occur, we fix p and find the associated critical temperature Tc{p) 
for each {L,L'). As usual, we take periodic boundary conditions along the finite direction 
of the strips, and L' = L — 1. Sparse-matrix techniques were used for the diagonalisation 
of the thermal transfer matrices. We have obtained the approximate phase diagrams using 
phenomenological scaling with L up to 7. For each temperature and concentration at the next 
step, L = 8, one would need to take into account nearly 32,000 configurations of connected 
pairs of columns and extract the two largest thermal eigenvalues of the respective transfer 
matrices; though symmetry considerations can reduce the number of distinct configurations 
involved, analysis of a few selected points convinced us that no additional insight would be 
gained this way. Data from L = 5, 6 and 7 scalings are shown in Fig. 0, together with an 
interpolated curve incorporating the exact values of Tc{p = 1) and of the limiting slope at 
p = 1, as well as the best estimate for the percolation threshold pc = 0.592745 ± 0.000002 



22 1 . The overall picture is similar to those provided by simple analytic scaling approaches 



11 and Monte Carlo simulations 25 



At the extreme points p = Pc, T^ = and p = 1, T^ = Tc{p = 1) it is known (respectively 



from Refs. |]T2| and [0]) that the approach of strip-scaling fixed points to the exact critical 
parameters is monotonic as L grows; as remarked in (a) and (b) above, these sequences are 
automatically reobtained in our scheme. For intermediate concentrations, the same uniform 
convergence does not seem to be present, as curves from successive scalings cross each other. 
This is almost certainly related to the approximations involved in the truncations referred 
to above. However, the curves are very close to each other, and to the interpolated phase 
boundary. Thus, although the truncation effects do not die away with increasing L (and 
there is no a priori reason they should do so), they are most probably of a quantitatively 
limited nature. 

We now turn to results for the reduced slope (1/Tc(l)) {dTc{p) / dp) ^^ of the critical curve 
at the pure limit, which is known exactly |T^. Close to p = 1 we were able to reach L = 11 



by using a truncated basis of connected states, consisting only of pairs of adjacent columns 
with at most two vacant sites overall. We checked the consistency of results thus obtained 
against those from the full basis (for L > 6) and from a basis with up to four vacant sites 
per two-column configuration (for L > 8). We made use of analytical expressions linking 
dTc{p)/dp with d^/dT and d^/dp, through the condition ^l/L = ^l'/L', which defines the 
critical line. Both partial derivatives can be evaluated by first-order perturbation theory 
instead of numerically; thus, calculations were performed very close to p = 1 {1 — p = 10~^ 
was typical), without loss of accuracy . 

Results from (L, L — 1) scalings up to L = 11 are given in Table 1, which complements the 
corresponding table of Ref. 0]. The sudden discontinuity when one goes from (L, L')= (7,6) 
to (8,7) is similar to what takes place for the fixed point of percolation between (L, L') = 
(4,3) and (5,4) [0,^], where convergence is interrupted. The extrapolation quoted in Table 
1 was obtained by taking the last four points of the sequence, and searching for the value of 



ip such that the plot of those data against L'"^ was the best possible straight line (see the 
work of Derrida and Stauffer in Ref. ||12|). Here, one has ip ~ 2.0. The exact slope is due to 
Thorpe and McGurn [|14| . The behaviour of the critical curve near the percolation threshold 



is determined by a crossover exponent giving the power-law dependence of e"^'^/^'"*^^^ on 
(p — Pc)- As (p is known exactly to be 1 [^, this can provide another test of the quality of 
the strip scaling results. 

At low temperatures, our ansatz is asymptotically exact. Care must be taken, however, 
with the spin degrees of freedom of occupied sites that are not directly connected to the 
origin |]T^,^|. With periodic boundary conditions, these may occur for L > 4. No matter 



whether they belong to finite clusters or are connected by a path that goes forwards before 
turning backwards to the origin, such sites do not contribute to the spread of thermal 
correlations along the backbone of the infinite cluster (they do, on the other hand, matter 
for the statistics of geometric configurations, so their probabilistic weight must be correctly 
included when counting the latter). In the corresponding thermal transfer matrices, the 
occupied unconnected sites are then treated as if they were absent. This is done, of course, 
for all temperatures and concentrations, but has well-marked effects in the low-temperature 
region. We have found that, if unconnected spins are not frozen out, the Tc x p curves 
for L > 4 develop an unphysical overhang ai T —>■ 0, p < Pc before homing in towards the 
correct Pc at T = 0. Once these corrections are incorporated, we get the correct vertical 
approach to Pc in the T — p plane. 

Slow convergence problems put a limit to the lowest temperatures within reach of inves- 
tigation; this is usually of order 10% of the critical temperature for the pure system, Tc(l); 
for the smallest strip widths, / = 2 - 4 we get to ~ 0.06Tc(l). As the phase boundary is 
vertical in that neighbourhood, the corresponding values oip—pc are ~ 10~^. This is enough 
to see the exponential behaviour, with = 1 along more than two orders of magnitude of 
variation in p — pc- In Fig. § we show, for L= 5, 6 and 7, the quantity 

_ exp(-2J/re(p)) , 

<--L = 7-r\ — ' v^J 

against ln(p — pc{L)) (where Pc{L) is the fixed point of phenomenological renormalisation, 
and Tc{p) is obtained from our ansatz, both for strips of widths L and L — 1). For L= 6 
and 7, the region in which Cl is approximately constant goes from p — Pc{L) ~ 4 x 10~^ to 
~ 6 X 10^'^. The approximate values at which the Cl become stable are: 6.1, 3.8 and 2.3 
respectively for L= 5, 6 and 7. It does not seem possible to obtain a reasonable estimate 
of the limiting value of this coefficient. Coo from extrapolation of this sequence. For bond 



dilution, a value for comparison is p4[ 21n2 = 1.386 .... A heuristic generalisation for site 
dilution (see |]1T|, pg. 190) would lead to Coo = (1/Pc)ln2 = 1.169.... The interpolated 
curve used in Fig.^ has C ~ 4, for consistency with corresponding data of Ref. [l^; for 
purposes of overall comparison of the shape of the phase diagram, the difference relative to 
the presumed value ~ 1.169 is not relevant. 

It must be remarked that the region covered in Fig.|^ is very difficult to reach e.g. in 
Monte-Carlo simulations, as one is at concentrations extremely close to criticality. For 



comparison, the lowest p used in recent work [25| corresponds to p — Pc — 7 x 10 ^. From 



the data depicted above, it seems safe to conclude that the exponential behaviour is dominant 
at least up to a critical temperature of order Tc/J ~ 0.7. 



Conformal invariance p6[ allows one to extract additional information from strip scaling, 
via the relationship between correlation- length amplitudes on a strip of width L at criticality 
and the decay-of-correlations exponent tj: 

V = L/naTc) . (7) 

In the present case, one has two questions to answer: first, whether conformal invariance still 
is valid for random systems: second, if so, how to define the correlation length which enters 
Eq.|^. As for the former, transfer-matrix and field-theoretical results [^] indicate the 
affirmative, provided that one considers averages over disorder. The latter is more involved, 
with the two most obvious candidates being the correlation decay factors related to the 



"most probable" and "average" correlation functions [^6|,|T7|. Based on numerical evidence 
for the bond-disordered two-dimensional Ising system |19|, we assume that for dilution the 
distinction between results coming from either definition amounts at most to logarithmic 
corrections. Thus, we shall use the average correlation length given by Eq.^ which, as 
remarked above, is closely related to an average of the two largest Lyapunov exponents, i.e., 
the most probable correlation decay. 

In Fig.^ we display the values of t] along the approximate critical curves. As can be 
seen, behaviour is not uniform with L; however, the overall trend for the exponent is to keep 
an approximately constant value close to that of the pure Ising model, rji = 1/4, and drop 
towards the percolation value rjp = 5/24 close to pc- 

The above-mentioned trend becomes much more apparent, and uniform with L, when 
one considers the variation of rj along the interpolated critical curve. Fig.^ illustrates this. 
Note the difference between vertical scales in Figs.^ and ^. Recall from Fig.Q that the 
approximate and interpolated curves are very close in T — p space. This shows that the 
averaged correlation length is very sensitive to slight variations in the parameters. In both 
figures, it is apparent that the values oi rj at p = pc, p = 1 approach the exact values 
as L increases. Again, the two special cases of the pure Ising model and percolation have 
been treated previously [T^ll^, and our results for these limits coincide with those already 



obtained. 

The behaviour depicted above is to be compared with current field theory results for 
random Ising models in two dimensions p7|-p9|. The theories agree in predicting that the 
specific heat singularity is of In ln(T — T^) type rather than the ln(T — T^) of the pure case. 
But one group of theories ^Tf predicts that the asymptotic pair correlation is as in the pure 



case (< o-Qa-r > ~ r"'', rj = 1/4) while the other p8[ concludes that the correlation has 
the form < aoaj. > ~ g-A{inr)2^ 

Our results clearly support the first class of predictions, and show no sign of the huge 
change of t] expected if the second prediction were correct. 

Large Monte Carlo simulations on the random-bond Ising model at a particular (self- 
dual) concentration p9| also provide indirect evidence in favour of the first class of results. 
An evaluation of the susceptibility exponent 7 is taken together with scaling relations to 
infer a value for rj, which is close to the pure value 1/4. A direct evaluation of r], again 
pointing to the pure result, has also been obtained for the same random bond Ising 
model by strip scaling on long realisations of the random system. Very recent Monte Carlo 



results for the correlation functions at criticality also support this view pn|. 

In conclusion, application of the calculational scheme just described to the site-diluted 



Ising model gives reliable and fairly accurate results, though convergence towards the exact 
values as strip width L increases seems to be generally non-uniform. Our data for the expo- 
nent 7] at intermediate concentrations give clear evidence in favour of one of the competing 



classes of field-theoretic treatments p8| , p7| , p9| , p0 



A straightforward extension of the method to the random-field Ising model in two di- 
mensions gives results for the averaged correlation length which are qualitatively similar to 
those predicted from the two largest Lyapunov exponents [^. These, in turn, seem to differ 
appreciably from the correlation lengths obtained directly from direct fits of the averaged 
correlation decay. As remarked above, cases such as this where frustrations are present tend 
to much subtler than when all correlations are ferromagnetic. We defer a detailed analysis 
of this point to a forthcoming publication. 
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TABLES 
TABLE I. Reduced slope at p = 1 



L/L' 


1 (dTAp)\ 
Tc(l) \ dp Jp^i 


3/2 


1.4461 


4/3 


1.4765 


5/4 


1.5032 


6/5 


1.5165 


7/6 


1.5223 


8/7 


1.4944 


9/8 


1.5015 


10/9 


1.5065 


11/10 


1.5101 


Extrapolated 


1.53 ± 0.015 


Exact 


1.565 



10 



FIGURES 

FIG. 1. Semi-logarithmic plot of correlation function < aoaji > (squares) against distance R 
for the pure Ising model on strips of width L of a square lattice. Values of ^ in lattice spacing 
units, are as given by Eq. 2. Straight lines have slope —1/S,- Tc is the critical temperature of the 
two-dimensional lattice. 

FIG. 2. Approximate phase diagram from (L, L — 1) scalings for L =5, 6 and 7 and from 
interpolation between exact limits. 



FIG. 3. Exponential behaviour of Tc against p close to pc- See Eq. 6 for definition of Cl- 

FIG. 4. Correlation exponent t] from correlation-length amplitudes, along approximate 
{L,L — 1) critical curves. Black squares are at (0.59275, 5/24) and (1, 1/4) . 

FIG. 5. Correlation exponent t] from correlation-length amplitudes, along interpolated critical 
curve. Black squares are at (0.59275, 5/24) and (1, 1/4). 
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